home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / fft / hc_init.c < prev    next >
Encoding:
Text File  |  2001-08-22  |  3.2 KB  |  129 lines
  1. /* fft/hc_init.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. TYPE(gsl_fft_halfcomplex_wavetable) *
  21. FUNCTION(gsl_fft_halfcomplex_wavetable,alloc) (size_t n)
  22. {
  23.   int status;
  24.   size_t i;
  25.   size_t n_factors;
  26.   size_t t, product, product_1, q;
  27.   double d_theta;
  28.  
  29.   TYPE(gsl_fft_halfcomplex_wavetable) * wavetable ;
  30.  
  31.   if (n == 0)
  32.     {
  33.       GSL_ERROR_VAL ("length n must be positive integer", GSL_EDOM, 0);
  34.     }
  35.  
  36.   wavetable = (TYPE(gsl_fft_halfcomplex_wavetable) *) 
  37.     malloc(sizeof(TYPE(gsl_fft_halfcomplex_wavetable)));
  38.  
  39.   if (wavetable == NULL)
  40.     {
  41.       GSL_ERROR_VAL ("failed to allocate struct", GSL_ENOMEM, 0);
  42.     }
  43.  
  44.   wavetable->trig = (TYPE(gsl_complex) *) malloc (n * sizeof (TYPE(gsl_complex)));
  45.  
  46.   if (wavetable->trig == NULL)
  47.     {
  48.       /* error in constructor, prevent memory leak */
  49.  
  50.       free(wavetable) ; 
  51.  
  52.       GSL_ERROR_VAL ("failed to allocate trigonometric lookup table", 
  53.             GSL_ENOMEM, 0);
  54.     }
  55.  
  56.   wavetable->n = n ;
  57.  
  58.   status = fft_halfcomplex_factorize (n, &n_factors, wavetable->factor);
  59.  
  60.   if (status)
  61.     {
  62.       /* error in constructor, prevent memory leak */
  63.  
  64.       free(wavetable->trig) ; 
  65.       free(wavetable) ; 
  66.  
  67.       GSL_ERROR_VAL ("factorization failed", GSL_EFACTOR, 0);
  68.     }
  69.  
  70.   wavetable->nf = n_factors;
  71.  
  72.   d_theta = 2.0 * M_PI / ((double) n);
  73.  
  74.   t = 0;
  75.   product = 1;
  76.   for (i = 0; i < n_factors; i++)
  77.     {
  78.       size_t j;
  79.       const size_t factor = wavetable->factor[i];
  80.       wavetable->twiddle[i] = wavetable->trig + t;
  81.       product_1 = product;    /* product_1 = p_(i-1) */
  82.       product *= factor;
  83.       q = n / product;
  84.  
  85.       for (j = 1; j < factor; j++)
  86.     {
  87.       size_t k;
  88.       size_t m = 0;
  89.       for (k = 1; k < (q + 1) / 2; k++)
  90.         {
  91.           double theta;
  92.           m = m + j * product_1;
  93.           m = m % n;
  94.           theta = d_theta * m;    /*  d_theta*j*k*product_1 */
  95.           GSL_REAL(wavetable->trig[t]) = cos (theta);
  96.           GSL_IMAG(wavetable->trig[t]) = sin (theta);
  97.  
  98.           t++;
  99.         }
  100.     }
  101.     }
  102.  
  103.   if (t > (n / 2))
  104.     {
  105.       /* error in constructor, prevent memory leak */
  106.  
  107.       free(wavetable->trig) ; 
  108.       free(wavetable) ; 
  109.  
  110.       GSL_ERROR_VAL ("overflowed trigonometric lookup table", GSL_ESANITY, 0);
  111.     }
  112.  
  113.   return wavetable;
  114. }
  115.  
  116.  
  117. void
  118. FUNCTION(gsl_fft_halfcomplex_wavetable,free) (TYPE(gsl_fft_halfcomplex_wavetable) * wavetable)
  119. {
  120.  
  121.   /* release trigonometric lookup tables */
  122.  
  123.   free (wavetable->trig);
  124.   wavetable->trig = NULL;
  125.  
  126.   free (wavetable);
  127. }
  128.  
  129.